Energetically stable particle-like Skyrmions in a trapped Bose-Einstein condensate 
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We numerically show that a topologically nontrivial 3D Skyrmion can be energetically stable in 
a trapped two-component atomic Bose-Einstein condensate, for the parameters of 87 Rb condensate 
experiments. The separate conservation of the two atomic species can stabilize the Skyrmion against 
shrinking to zero size, while drift of the Skyrmion due to the trap-induced density gradient can be 
prevented by rotation or by a laser potential. 
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Localized topological excitations that do not perturb 
the order parameter field at large distances from the par- 
ticle, and which are characterized by a topologically in- 
variant winding number, are well-known in nuclear and 
elementary particle physics 0,0 While their study 

in nuclear physics remains an experimental challenge, 
the recent experimental progress in atomic Bose-Einstein 
condensates (BECs) with internal spin degrees of freedom 
0, Q has raised the possibility of the existence of well- 
localized topological Skyrmions in atomic gases 0, @> • 
In this Letter we identify, and show how to overcome, the 
specific instabilities of Skyrmions in trapped two-species 
atomic BECs, and hence demonstrate their energetic sta- 
bility under realistic experimental conditions. 

Battye et at [t| recently considered an infinite homoge- 
neous two-species BEC, with constant total atom density. 
They showed that an energetically stable Skyrmion may 
exist as a result of phase separation of the two species, 
which suppresses the decay. These calculations were ex- 
tended to non-constant total atom density, and to the 
trapping of one component [loj . In this paper we show 
that in a harmonically trapped system there are addi- 
tional instabilities, not considered in Refs. 0,0], which 
will play a crucial role in the experimental realization of 
Skyrmions in atomic BECs. Additional physical mech- 
anisms, such as rotation or optical potentials, will be 
required for stability. We also show how densit y fl uctu- 
ations, associated, e.g., with phonon emission are 
important in the Skyrmion decay process. 

There has been an explosion of interest in vortex and 
soliton experiments in atomic BECs |l2|. and we antici- 
pate similar developments for other topological objects. 
Hence we identify the Skyrmion energetic stability cri- 
teria for the parameters of the JILA two-species 87 Rb 
experiments |SJ. We find a threshold frequency below 
O.lw, when only one species is rotated, and a narrow win- 
dow of rotation frequencies for the entire system around 
0.085cj. We numerically evaluate the stable configura- 
tions (Fig. ^) by minimizing the energy of the full 3D 
mean-field theory of coupled Gross-Pitaevskii equations 
(GPEs) and calculate the associated topological charges. 



A Skyrmion is a topological particle-like soliton solu- 
tion with a coreless 3D texture. Besides their intrinsic 
fundamental interest, Skyrmions have important appli- 
cations in nuclear physics 0, 0, Q , and analogous struc- 
tures are postulated for early Universe cosmology [l3| . 
Skyrmions are localized objects such that the order pa- 
rameter field becomes uniform sufficiently far from the 
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FIG. 1: Density and order parameter profiles for energetically 
stable trapped Skyrmions (N+ = N- = 4.5 x 10 e ). Top: 3D 
densities. The central (blue) torii are isosurfaces of \ip- \ 2 . Iso- 
surfaces of |V>+| 2 (red), are shown for x < 0: on the y-z plane 
between the isosurface sections its density is indicated by a 
colormap from red (lowest) to purple (highest). Left: Stabi- 
lized by rotating ip- only, with angular velocity 0.1a;. Right: 
Stabilized by rotating the entire system with angular velocity 
0.085oj. Bottom: ID densities (left axis, units of x^) and 
the order parameter profile X(x, 0, 0) (right axis, dotted line) 
for rotating ip- only. Solid line: \ip+\ 2 + \4>-\ 2 ■ Dashed line: 
ji/i+j 2 . Dash-dot line: I'i/'-l 2 - A(a;,0, 0) is extracted from the 
wave functions. The notch in A is due to the mean-field repul- 
sion inflating the vortex ring core. A(r) is highly anisotropic, 
with A(0, 0, z) qualitatively close to 2 arctan(|z|). 
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particle. As a result, we may enclose the Skyrmion by 
a sphere on which the field configuration has a constant 
asymptotic value. This sphere can be mapped to a single 
point in the order parameter space. Topologically, it indi- 
cates that the 3D physical region inside the sphere can be 
represented by a compact 3-sphere S 3 (a sphere in 4D) 
UJ. Here we consider a SU(2) order parameter space 
which topologically also corresponds to S 3 . Then the 
mappings from the compactified 3D physical region to the 
order parameter space are represented by the field U(r): 
S 3 — > SU(2). The crucial point is that having associated 
the enclosing sphere with a single point in S 3 , we place 
the elements of the physical space into a correspondence 
with the elements of the compact order parameter space 
SU(2). The mappings from S 3 to SU(2) fall into topo- 
logical classes, each characterized by an integer- valued, 
topologically invariant winding number W |lj, |3, li| : 

W = ^J d 3 rTr[U(d a U^U(d^)U(d u U^)} . (1) 

Here e a f}v denotes the completely antisymmetric tensor, 
the repeated indices are summed over, and di represents 
the derivative with respect to the spatial coordinate x%. 
The topological charge describes how many times the set 
of points in physical space is 'wrapped' over the order 
parameter space of SU(2) for a given element U(r). 

With cold atoms an SU(2) order parameter space is 
afforded most simply by a two-component BEC, whose 
interactions effectively fix the local value of the total den- 
sity |'0+(r)| 2 + |V'-( r )| 2 = p(r) of the two complex macro- 
scopic wave functions ip+ and 7\. We may then ex- 
press the two-component BEC as 

(^O-VWJIT'WQ. (2) 

We search for a topologically nontrivial solution for U (r) 
with a nonvanishing winding number, determined by 
Eq. Q . By minimizing the energy of the corresponding 
trapped BECs, we may investigate the energetic stability 
of the Skyrmion and find its equilibrium configuration. 
As an initial state for the numerical simulations we use 
U(t) = exp[iA(r)<7 • f] [J^, or a similar state displaced 
from the trap center, where o~i denote the Pauli spin ma- 
trices and f represents the unit radial vector. A direct 
substitution into Eq. (TJ yields W — 1 for a monotonic 
function A with A(0) = and A = n at the gas boundary. 
The corresponding BEC wave functions read: 

AM r )\ _ rrfcos[X(r)} -isin[A(r)]cos0\ , , 
\ip-(r)J -vPy -isin[A(r)]siii0e X p(#) J ' { > 

Here (A, 8, <j)) can also be understood as the spherical 
angles of the 3-sphere. Then the boundary of the atom 
cloud corresponds to the pole A = n of the 3-sphere. 

Due to the topological stability of the Skyrmion, any 
continuous deformation of Eq. @, without altering the 



asymptotic boundary values, is still a Skyrmion with 
W = 1 . Here, ip- (the "line component" ) forms a vortex 
line with one unit of circulation around the core oriented 
along the z axis (Fig. G). Moreover, ip + (the "ring com- 
ponent" ) forms a vortex ring with one unit of circulation 
around its core on the circle z — 0, r — A _1 (^) [3. The 
line component is confined in a toroidal region inside the 
ring component in the neighborhood of the ring core. 

Despite its topological stability, the Skyrmion solu- 
tion in several physical systems is energetically unstable 
against shrinking to zero size without additional stabiliz- 
ing features. This can be understood by means of simple 
scaling arguments: if the kinetic energy density, or the 
order parameter 'bending energy', is quadratic in the gra- 
dient of the order parameter, this energy density scales 
as 1/R 2 with respect to the size R of the Skyrmion. Be- 
cause the Skyrmion occupies a volume proportional to 
R 3 , the energy E cx R monotonically decreases with the 
size of the Skyrmion. In the Skyrme model, the stability 
is provided by an additional interaction term, with the 
energy density scaling as 1/i? 4 0. 

It was proposed that, in a two-component BEC, the 
separate conservation of the species can effectively sta- 
bilize the Skyrmion in a homogeneous space against col- 
lapse to zero size . In the regime of phase separation, 
with the scattering lengths satisfying a ++ a__ < a 2 + _ 1 
the two species can strongly repel each other and the 
toroidal filling due to the line component can prevent the 
vortex ring from shrinking to zero radius 0>lll- Conse- 
quently, a filled vortex ring, as opposed to an empty vor- 
tex ring, can be more stable against collapse. Moreover, 
in the Skyrmion the filling has one unit of angular mo- 
mentum about the z axis resulting in a 1/r 2 centrifugal 
barrier that further prevents the shrinking. This should 
be contrasted to the Skyrmions in a ferromagnetic spin- 1 
BEC 0, which are closely related to the Skyrmion ex- 
citations proposed in superfluid liquid 3 He |f6j, and can 
freely mix the atoms between the different spin compo- 
nents. As a result, no stabilizing features due to the atom 
number conservation exist in that case. 

We next turn to numerical studies of the energetic sta- 
bility of the Skyrmion using the full 3D mean-field theory 
of the coupled GPEs without additional simplifications: 

itnfri = ( - |^V 2 + V(r) + £ K ife |^| 2 )^ • (4) 

k 

Our simulations are performed for the parameters of 
the JILA two-species vortex experiment using perfectly 
overlapping isotropic trapping potentials for both states 
V(r) — muj 2 r 2 /2, with u> — 2ir x 7.8 Hz and harmonic 
oscillator length xi lo = (h/muj) 1 / 2 ~ 3.86//m In the 
interaction coefficients, = ^■nh 2 aijNi/m, the number 
of atoms in each species is represented by N%, an denotes 
the intraspecies, and (i 7^ j) the interspecies, scatter- 
ing length. For the \F — 2, mj = 2) and |1, f) spin states 
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of 87 Rb we have [l7|: a ++ = 5.67nm, a = 5.34nm, 

= 5.50nm. The GPEs depend on the dimensionless 

interaction coefficients n'ij = AirNidij jx\- lo . Hence, the 
dynamics is unchanged in any scaling of length and time, 
which does not change Nfu>. 

We found the ground states by imaginary time evolu- 
tion of the GPEs igj, using Runge-Kutta 18] and split- 
step algorithms. At every time step, we separately nor- 
malized both wave functions to fix the atom number in 
each component. The numerical simulations were fully 
3D, allowing for cylindrically asymmetric dynamics. The 
most demanding numerics was performed on a parallel 
multiprocessor supercomputer, using up to 32 processors. 
Spatial grids of 128 3 , 256 3 , and 512 3 were used. For 
a typical spatial range of about 30x^ o the correspond- 
ing grid spacings are respectively 0.23x^ o , 0A2xh o , and 
O.OGxho- A typical imaginary time step was 0.0025/w. 

As a cylindrically symmetric initial state we used the 
Skyrmion (J3J) with the Thomas- Fermi density profile , 
resulting in the winding number W = 1. The winding 
number was calculated during the imaginary time evo- 
lution to determine the stability of the Skyrmion. For 
small BECs, with N T = 10 4 atoms (N T = N+ + N-), 
the nonlinear repulsion between the two species was too 
weak to prevent the Skyrmion from shrinking for any 
N-JNt- Also the line component ^_ rapidly diffused 
to the boundary of the finite-size atomic cloud, altering 
the boundary conditions and resulting in a decreasing 
winding number W < 1. Because the topological sta- 
bility of the Skyrmion necessitates a well-defined asymp- 
totic boundary condition with no phase variation, the 
Skyrmion was clearly lost. This decay mechanism is char- 
acteristic of trapped BECs and cannot occur in homoge- 
neous systems, which implicitly assume large N + /Nt- 

For large BECs with Nt > 10 6 , represented in Fig.|2 
the same instability mechanism was observed with a very 
large fraction of line component N^/Nt > 0.75. In other 
cases displayed in Fig. the winding number typically 
evaluated to either or 1 to better than 1%, with A = 7r 
well preserved at the boundary. A drop from W = 1 
to W = indicated instability of the Skyrmion against 
shrinking, see inset to Fig. EI The Skyrmion collapse via 
shrinking occurred when the high density central part of 
the ring component, which passes through the line com- 
ponent vortex core, pinched off. The total density varies 
significantly at the trap center when the vortex ring core 
collapses, emphasizing the importance of density fluctu- 
ations in the decay process, see Fig. El (top). With suf- 
ficiently large BECs, even for the collapsed Skyrmions 
with no vortex ring, the local energetic minimum can 
correspond to line component trapped inside a toroidal 
region with ip + forming an atom-optical confinement. 

For sufficiently large total and line component atom 
numbers, the nonlinear repulsion between the two species 
was strong enough to inhibit the collapse of the vortex 
ring, stabilizing the Skyrmion against shrinking for cylin- 
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FIG. 2: Skymion stability diagram. Total atom number Nt 
vs. the fraction in the line component N-/Nt- X indicates 
instability. © indicates stability of the cylindrically symmetric 
state against shrinkage, which may be stabilised against drift 
as discussed in the text. 



drically symmetric initial states (Fig. However, due 
to the inhomogeneous potential, the Skyrmion was still 
unstable with respect to drift towards the edge of the 
BEC where its energy is lower. This was found using 
cylindrically asymmetric initial states, emphasizing the 
importance of the full 3D simulation. The drift occurred 
with the vortex line moving towards the boundary. Once 
the vortex line drifted to a low-density region, the nonlin- 
ear repulsion was no longer strong enough to prevent the 
shrinking of the Skyrmion and the collapse occurred as 
described above, see Fig. |3| (bottom) . The drift reduces 
the total angular momentum of the atoms indicating an 
energetic instability. Physically, this results from dissi- 
pation, e.g., due to thermal atoms. 

We have investigated several mechanisms by which 
Skyrmions might be stabilized against the vortex drift 
instability. Perhaps the simplest is to rotate the line vor- 
tex component about the z axis, above its critical angular 
velocity, thus stabilizing the vortex. This should always 
work for cylindrically symmetric traps in which the sym- 
metric density prevents the rotation coupling into the 
other BEC component. For the case of Fig.^ a hne com- 
ponent angular velocity of O.Ilu stabilised the Skyrmion. 

Rotation of both species introduced a new instability 
mechanism. For sufficiently rapid rotations the line com- 
ponent again reached the boundary, resulting in W < 1. 
At the same time the ring vortex accomodated the rota- 
tion by twisting and finally breaking at the gas bound- 
ary. Nevertheless, for the case of Fig.^ we found a small 
range of angular velocities around 0.085u; for which the 
Skyrmion was stable against this, and against line vor- 
tex drift. It stabilized off the rotation axis, indicating an 
equilibrium between the lower energy of the line vortex 
and the higher energy of the ring component threading 
through it with increased distance from the trap center. 
The breaking of cylindrical symmetry implies a family of 
degenerate Skyrmions parametrized by the polar angle <j>. 

Another method for stabilizing the Skyrmion is to in- 
hibit the drift towards low-density regions by creating a 
positive density gradient around the vortex line. This 
can be implemented with a blue-detuned Gaussian laser 
beam along the z axis |l9| , providing a cylindrically sym- 
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called Skyrmions |2l|. Although these vortices may be 
energetically stable in a rotating trap j2^| , they are topo- 
30 logically trivial in the 3D ferromagnetic spin-1 BEC and, 
analogously to the vortex lines with two units of circu- 
lation, they can be continuously deformed to a topologi- 
cally trivial uniform spin distribution. 

This research was supported by the EPSRC, the Aus- 
tralian Research Council, and the Australian Partnership 
for Advanced Computing. 



FIG. 3: Decaying Skyrmion's imaginary time evolution. Top: 
Densities (units of x^f) along the x axis. |V>+| 2 : it = 20 
(dash-dot), 27.75 (dashed), 28 (solid), 50 (dotted). The dot- 
ted curve which is zero at x = is 1 2 at it — 50: for other 
times it is qualitatively similar. W — by it = 28.25, and 
it = 50 is the stationary state. N+ = 7V_ = 2.3 x 10 6 . Inset: 
Winding number W versus imaginary time. Bottom: Densi- 
ties on the z — plane, left IV 1 -! 2 ! right |^/>+| 2 . Since there 
is no stabilisation, the line vortex core (red circle) has moved 
towards the boundary and the ring vortex (surrounding the 
blue circle) is about to collapse. Colormap and parameters as 
for Fig. 1. 

metric repulsive Gaussian dipole potential perpendicular 
to the z axis. The total potential then has the "Mexican 
hat" form: V = muj 2 r 2 /2 + V exp[-2(x 2 + y 2 )/w 2 ]. Set- 
ting w = 7xho and Vb = 25hio, corresponding to about 
twice the width of the vortex line core, successfully sta- 
bilized a Skyrmion with N + = N_ = 8 x 10 6 . 

Proving numerically that a Skyrmion is stable is diffi- 
cult. The essential requirement is that it be stationary 
under imaginary time evolution. This was determined by 
plotting the density on ID sections through the system, 
and the phase variation on 2D slice planes. In unsta- 
ble cases these would evolve until the Skyrmion decayed. 
When they became asymptotically constant as a function 
of imaginary time the Skyrmion was judged to be stable. 
Convergence was particularly slow, and hence difficult to 
establish, when both components were rotated. 

The Skyrmion can be created using electromagnetic 
fields to imprint topological phase singularities on the 
matter field while changing the internal state of the 
atoms, as proposed in Ref. 7J. In particular, a vortex 
ring can be engineered in a controlled way with an ap- 
propriate phase-coherent superposition of three orthog- 
onal standing waves. Due to dissipation, the prepared 
Skyrmions relax to the ground states calculated here. 

Different spin textures are often referred to as 
"Skyrmions" , even when they are not characterized by 
the topological invariant (JIJ. It is important to em- 
phasize that the Skyrmions studied in this paper are 
fundamentally different from the nonsingular Anderson- 
Toulouse or Mermin-Ho vortices [2(J, frequently also 
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